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Abstract 

The  purpose  of  the  algorithm  developed  in  this  thesis  was  to  create  a  post  processing 
method  that  could  resolve  objects  at  low  signal  levels  using  polarization  diversity  and  no 
knowledge  of  the  atmospheric  seeing  conditions.  The  process  uses  a  two-channel  system, 
one  unpolarized  image  and  one  linearly  polarized  image,  in  a  GEM  algorithm  to 
reconstruct  the  object.  Previous  work  done  by  Strong  showed  that  a  two-channel  system 
using  polarization  diversity  on  short  exposure  imagery  could  produce  images  up  to  twice 
the  diffraction  limit.  In  this  research,  long  exposure  images  were  simulated  and  a  simple 
Kolmogorov  model  used.  This  allowed  for  the  atmosphere  to  be  characterized  by  single 
parameter,  the  Fried  Parameter.  Introducing  a  novel  polarization  prior  that  restricts  the 
polarization  parameter,  it  was  possible  to  determine  the  Fried  Parameter  to  within  half  a 
centimeter  without  any  addition  knowledge  or  processes.  It  was  also  found  that  when 
high  polarization  diversity  was  present  in  the  image  could  be  reconstructed  with 
significantly  better  resolution  and  signal  level  did  not  affect  this  resolving  capability.  At 
very  low  signal  levels,  imagery  with  low  to  no  diversity  could  not  be  resolved  at  all 
whereas  high  diversity  resolved  equally  as  well  as  if  there  was  a  high  signal  level. 

Current  algorithms  being  used  do  not  include  polarization  diversity  but  can  substantially 
improve  resolution.  Application  of  this  algorithm  could  be  used  in  dim-object  detection 
around  satellites  as  well  as  solar  surface  imagery. 
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BLIND  DECONVOLUTION  THROUGH  POLARIZATION  DIVERSITY  OF 

LONG  EXPOSURE  IMAGERY 


I.  Introduction 


1.1  Motivation 

In  200 1  the  Department  of  Defense  released  a  comprehensive  report  on  the  United  States 
Space  Capabilities.  In  that  report,  it  was  said  that  we  are  ripe  for  a  “Space  Pearl  Harbor.” 
[3]  Since  then,  there  has  been  a  concerted  effort  to  mitigate  this  possibly  with  the 
advancement  of  Space  Superiority.  This  is  broken  down  into  three  categories:  Offensive 
Counterspace,  Defensive  Counterspace,  and  Space  Situational  Awareness  (SSA).  In  orbit 
around  the  earth  it  is  very  difficult  to  identify  and  characterize  anomalies  that  may  occur 
with  “blue”  spacecraft  or  the  functions  and  purpose  of  “red”  spacecraft.  That  is  where 
SSA  comes  in.  It  is  the  attempt  to  have  complete  awareness  of  the  battlespace  in  orbit. 

Advanced  sensors  designed  to  inspect  the  orbital  battlespace  or  ground-based  telescope 
systems  are  required.  The  design  and  launch  of  satellites  are  very  costly,  especially  at 
geosynchronous  (GEO)  orbit.  At  Geo,  there  is  so  much  distance  between  satellites  that 
space-based  optical  systems  need  to  be  maneuvered  close  to  each  Resident  Space  Object 
(RSO)  of  interest.  This  greatly  limits  lifetime  due  to  finite  fuel,  and  also  restricts  the 
response  time  kill  chain  after  an  event  occurs.  On  the  other  hand,  ground  based  optical 
systems  have  a  comparably  low  cost,  can  be  easily  repaired  or  upgraded,  and  can  respond 
quickly  when  an  event  occurs.  The  drawback  is  that  observing  must  take  place  anywhere 
from  hundreds  of  kilometers,  for  low  earth,  to  thousands  of  kilometers  for  GEO.  On  top 
of  that,  resolution  is  reduced  considerable  by  atmospheric  seeing  conditions.  If  large 
enough  telescopes  or  telescope  arrays  are  constructed,  atmospheric  distortion  is  the  main 
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thing  that  needs  to  be  mitigated.  Adaptive  Optics  can  help  significantly  with  this  but  do 
not  correct  for  all  atmospheric  distortion.  Post  processing  algorithms  can  be  used  to 
further  reconstruct  the  RSO.  Combinations  of  Adaptive  Optics  and  post  processing  are 
currently  being  used  operationally  to  characterize  satellites  and  anomalies  in  space. 


1.2  Goals 

1.  The  purpose  of  this  thesis  research  is  to  develop  an  algorithm  that  can  be  used  with  or 
without  adaptive  optics  to  improve  image  resolution  of  space  objects  through  the  use  of 
polarization  diversity. 

2.  The  algorithm  should  produce  better  resolved  images  when  polarization  diversity  is 
high. 

3.  Knowledge  of  the  atmospheric  seeing  parameter  should  not  be  necessary  to  restore  the 
image  and  further  should  be  capable  of  being  estimated  from  a  likelihood  model. 

4.  The  algorithm  should  be  able  to  function  even  at  very  low  signal  levels. 


1.3  Previous  Work 

The  worked  developed  in  this  thesis  is  built  primarily  from  the  research  done  by  Major 
David  Strong  [11]  and  Lieutenant  Colonel  Adam  McDonald  [7].  Strong’s  dissertation 
created  a  two-channel,  one  unpolarized  and  one  linearly  polarized,  blind  deconvolution 
algorithm  for  passively  illuminated  objects  as  is  done  in  this  thesis  but  with  several  key 
differences.  His  algorithm  was  created  for  use  in  short  exposure  images  as  opposed  to 
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long  exposure.  Details  of  the  object  will  not  be  blurred  out  as  much  from  averaging  over 
a  longer  period  of  time.  The  drawback  to  the  short  exposure  case  is  that  signal  levels  are 
significantly  reduced  and  therefore  the  SNR  is  much  lower.  The  other  downside  to  this  is 
that  the  point  spread  function  (psf)  of  the  atmosphere  is  not  as  well  known.  As  an  image 
is  integrated,  the  psf  will  tend  toward  a  well  known  and  easily  modeled  transfer  function, 
such  as  the  Kolmogorov  spectrum.  In  contrast,  if  the  integration  time  is  small, 
fluctuations  in  the  atmosphere  can  cause  the  psf  to  vary  greatly.  This  makes  it  difficult  to 
estimate  and  characterize  [11]. 

Another  major  difference  is  in  the  development  of  the  two-channel  algorithm  derivation. 
In  this  thesis,  a  prior  density  for  the  polarization  tenn  is  included  to  restrict  the  possible 
values  that  the  polarization  parameter  can  take.  This  allows  for  the  polarization  state  and 
the  object  to  be  estimated  simultaneously.  The  addition  of  the  prior  gives  significantly 
increased  information  when  polarization  diversity  is  present  and  when  calculating  the 
correct  seeing  parameter  of  the  atmosphere  [11  :Chap  5]. 

In  MacDonald’s  work,  he  developed  a  method  for  estimating  the  seeing  parameter  of  the 
atmosphere  (r0)  for  a  single  channel  system  of  laser  light,  conforming  to  a  negative 
binomial  distribution.  The  method  involved  representing  the  ro  probability  density 
function  in  some  distribution.  In  that  case,  based  on  the  observation  that  good  seeing 
(high  r0)  is  much  less  likely  than  bad  seeing  (low  r0),  an  exponentially  decreasing 
function  was  chosen 


/Oo)  = 


N2 

ravg 


e  r°-v9 


(1.1) 
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where  N  is  the  number  of  pixels  on  a  side  and  ravg  is  some  parameter  that  must  be 
iteratively  calculated.  Once  the  average  seeing  parameter  is  found  it  produces  a  graph 
similar  to  that  in  Figure  3.2.  When  equation  (1.1)  is  left  out  when  deriving  the  likelihood, 
there  is  no  decrease  around  the  actual  ro  value  and  the  likelihood  continues  to  increase 
forever  [7: Sect  3.2]. 

It  was  assumed  that  MacDonald’s  method  would  be  needed  to  find  the  correct  seeing 
parameter  for  the  algorithm  described  in  this  document.  Surprisingly,  this  was  not  the 
case.  By  using  the  polarization  diversity  algorithm  with  a  polarization  prior,  the 
likelihood  curve  naturally  produces  a  maximum  near  the  actual  ro  value.  The  curve  looks 
very  similar  to  what  is  produced  by  Adam  MacDonald’s  iterative  method. 


1.4  Summary  of  the  Document 

Chapter  1:  This  chapter  presents  a  brief  description  of  what  the  problem  with  current 
Space  Situational  Awareness  capabilities  and  why  it  is  important  that  new  methods  be 
developed  for  the  Air  Force.  Along  with  this,  relevant  previous  work  done  in  blind 
deconvolution  and  polarization  diversity  is  discussed. 

Chapter  2:  This  chapter  presents  the  background  material  necessary  to  understand  the 
development  of  the  polarization  diversity  blind  deconvolution  algorithm.  Topics 
discussed  include  polarization,  linear  systems,  Fourier  optics,  atmospheric  turbulence, 
and  estimation  theory. 

Chapter  3:  Development  of  the  blind  deconvolution  algorithm  is  explained  here. 
Beginning  with  a  two-channel  system,  Poisson  statistics  are  used  to  create  a  likelihood 
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model.  It  is  further  refined  by  introducing  a  complete  data  model.  The  general  method 
used  to  estimate  the  object  is  the  General  Expectation  Maximization  (GEM)  algorithm. 

Chapter  4:  Simulations  done  using  the  algorithm  and  its  results  are  discussed.  To 
quantify  how  the  algorithm  performs  at  varying  signal  levels  and  polarizations,  two  bar 
targets,  polarized  differently  in  the  same  scene,  are  propagated  through  a  simulated 
atmosphere  and  aperture  and  then  reconstructed. 

Chapter  5:  This  last  chapter  discusses  general  conclusions  about  the  functionality  and 
utility  of  the  developed  algorithm.  Also,  further  testing  and  application  are  mentioned  for 
future  work. 
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II.  Background 


2.1  Polarimetry 

2.1.1  Linearly  Polarized  Light 

Light  is  an  oscillating  electromagnetic  wave  with  both  an  electric  field  component  and  a 
magnetic  field  component.  The  entire  light  wave  can  be  described  completely  by  only 
the  electric  field  oscillations.  Given  a  cartesian  coordinate  plane,  the  total  electric  field  at 
any  point  is  given  by 

E(z,  t)  =  Ex (z,  t)  +  Ey(z,t )  (2.1) 

With  z  being  the  direction  of  propagation  and  t  is  time.  This  says  that  the  composite 
electric  field  is  just  a  superposition  of  the  x  and  y  components  of  the  field.  The 
individual  component  waves  are  described  by  the  well  known  standing  wave  equation 

Ex{z,  t )  =  Ex qCos  ( kz  —  (2-2) 

Ey(z,  t )  =  Ey0cos  (kz  —  a)t  +  e)j  (2.3) 

where  k  is  the  wave  number,  a)  is  the  frequency,  e  is  a  phase  shift  between  the  waves  and 
the  scalar  E0  is  the  magnitude  of  the  field  component.  Assuming  that  e  is  constant,  then 
the  total  electric  field  will  oscillate  linearly  at  some  angle  defined  by  the  phase  shift.  An 
example  of  this  is  shown  in  figure  2.1.  In  the  case  where  the  phase  shift  is  a  multiple  of 
n,  the  oscillation  will  be  along  either  the  x  or  y  axis.  This  is  called  horizontally  and 
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V 


Figure  2. 1 :  The  above  diagram  shows  a  total  oscillating  electric  field  that  is  vertically  polarized.  Image  taken  from  George  State 
University  Hyperphysics  website:  http://hyperphysics.phy-astr.gsu.edU/HBASE/phyopt/polclas.html#cl 


vertically  polarized  light,  respectively.  In  the  development  of  this  thesis,  only  linearly 
polarized  light  is  considered  so  elliptical  polarization  will  not  be  discussed  [9:280-285]. 


2.1.2  Unpolarized  Light  Sources 

The  term  “natural  light”  refers  to  a  source  that  is  completely  unpolarized.  However,  it  is 
not  that  the  light  is  unpolarized,  it  is  that  the  light  is  changing  its  polarization  state  so 
quickly  that  it  cannot  be  determined.  When  excited  electrons  in  an  atom  return  to  a  lower 
energy  state,  that  energy  is  released  via  a  photon  of  light.  The  electric  held  orientation  of 
the  emitted  photon  can  be  thought  of  as  random,  even  though  it  is  determined  by  the 
angular  momentum  of  the  system.  Photon  fluxes  can  vary  greatly  but  are  usually  on  the 
order  of  1020s-1.  With  such  a  large  number  of  constantly  changing  polarizations, 
natural  light  can  easily  be  considered  random  [9:280-285]. 
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2.1.3  Polarization  State 


For  this  thesis,  only  one  polarization  parameter  is  needed.  That  is  the  total  linear 
polarization  state,  P,  of  the  observed  light.  As  discussed  in  the  previous  section,  natural 
light  is  unpolarized.  When  that  light  is  integrated  over  time  and  decomposed  into  its  x 
and  y  components,  one  would  expect  to  find  half  of  the  light  to  be  horizontally  polarized 
and  the  other  half  to  be  vertically  polarized.  The  ratio  of  one  of  those  components  to  the 
total  light  collected  is  P.  It  is  unimportant  which  component  since  they  are  compliments 
of  each  other  and  coordinate  system  can  be  rotated  arbitrarily  in  the  x-y  plane.  In  order  to 
use  the  polarization  diversity  algorithm,  one  channel  collects  light  directly  off  the 
telescope  and  the  other  channel  utilizes  a  linear  polarizer  in  front  of  the  CCD.  The 
orientation  of  the  polarizer  does  not  matter  since  the  x-y  coordinate  system  chosen  is 
completely  arbitrary. 

When  natural  light  is  reflected  off  an  object,  it  can  become  polarized  by  the  material  the 
object  is  composed  of.  The  amount  of  polarization  is  dependent  on  material 
characteristics  and  the  geometry  of  the  scene.  Mamnade  objects  tend  to  polarize  light 
upon  reflection  very  strongly  as  do  sharp  edges. 


2.2  Diffraction 

2.2.1  Linear  System  Theory  &  The  Transfer  Function 

Linear  systems  are  those  that  produce  a  response  to  a  collection  of  inputs  equal  to  the 
sum  of  the  responses  obtained  from  the  inputs  if  they  were  applied  individually.  These 
types  of  transformations  can  be  found  in  many  places,  most  notably  in  the  time  domain 
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for  electronic  circuits  and  in  the  spatial  domain  for  optical  systems.  A  transfonnation 
operator,  L,  is  considered  to  be  linear  if  it  obeys  the  following  mathematical  property, 
defined  as  superposition 


L{ap(x)  +  bq(x )}  =  aL{p(x)}  +  bL{q(x)} 


(2.4) 


where  a  &  b  are  scalar  quantities  and  p{x)  &  q(x)  are  input  functions,  x  represents  an 
arbitrary  coordinate  system,  be  it  time,  space,  frequency,  etc. . .  This  superposition 
definition  has  some  important  properties  to  note.  First,  a  linear  system  is  unaffected  by 
the  scalar  magnitude  of  a  function  it  is  operating  on.  Secondly,  the  sum  of  the  output 
functions  is  equal  to  the  output  of  the  sum  of  the  functions;  meaning  that  there  is  no 
“mixing”  between  functions  when  the  operator  acts  on  several  at  one  time.  To  see  how 
this  applies  to  an  optical  system,  the  quantities  of  equation  (2.4)  must  be  given  physical 
definition  [6:19]. 


The  input  function  is  a  function  of  the  light  that  leaves  the  object.  It  has  not  passed 
through  any  optics  or  been  aberrated  in  any  way.  The  output  function  is  the  image 
obtained  by  the  camera.  The  linear  process  that  is  applied  to  the  object  is  the  actual 
propagation  of  the  light  through  the  atmosphere  and  optical  system.  To  understand  the 
fonn  of  the  linear  transformation,  the  impulse  response  of  the  system  must  be  found.  It  is 
assumed  known  that  if  a  linear  operator  is  applied  to  an  impulse,  the  output  is  the 
characteristic  response  of  the  system.  Mathematically,  an  impulse  can  be  described  by  a 
delta  function,  8  (x) .  Application  of  the  sifting  property  yields 


Pi  (*2) 


IXJ 

II 


p1(x1)8(x2  -x1)dx1 


(2.5) 
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Again,  x1  and  x2  are  two-dimensional  coordinates  and  p:  is  object  function.  If  the 
operator,  L,  is  applied  to  this  with  respect  to  the  x2  coordinate  space,  one  gets 


P2O2)  =  MPito)}  = 


II 


p1(x1)L{d(x2  -  x1)}dx1 


(2.6) 


As  mentioned  above,  the  impulse  response  is  defined  as  the  output  of  a  linear  operator  as 
applied  to  a  single  impulse.  This  will  be  given  the  definition 

h(x 2  -  xr)  =  L{d(x2  -  Xi)}  (2.7) 


substituting  this  back  into  equation  (2.6)  obtains 


P2(*2) 


00 

jj  Pi(Xi)/l(x2 
—  00 


—  x1)dx1 


(2.8) 


This  is  a  general  convolution  integral  between  the  original  object  function  and  the 
impulse  response.  Convolutions  are  almost  always  computationally  intense  and  difficult 
to  calculate.  If  this  is  done  in  the  spatial  frequency  domain,  the  two  functions  can  simply 
be  multiplied  together  and  then  the  inverse  Fourier  transform  taken  such  that 

p2(x2)=T~1{P1(fXi)H(fXi)}  (2.9) 
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Pi(/Xl)  and  H(fXi )  are  the  Fourier  transforms  of  pq  )  and  h(x2  —  xj,  respectively. 

is  the  inverse  Fourier  transform  operator.  Specifically,  H(fXi)  is  called  the 
Transfer  Function  of  the  system  and  is  defined  as  [6:19-21] 


00 

H(fX,fy )  =  //  h(x,y)e  l2n(fxX+fyy) dxdy 
—  00 


(2.10) 


2.2.2  The  Point  Spread  Function 

A  delta  function  was  used  to  define  the  impulse  response  of  the  system.  In  imaging 
applications,  the  impulse  response  is  also  referred  to  as  the  Point  Spread  Function  (psf). 
In  tenns  of  an  optical  system,  the  best  way  to  model  the  psf  is  to  use  a  completely  black, 
zero  intensity  image  with  one  pixel  at  the  center  having  a  maximum  intensity.  It  is  not  a 
perfect  delta  function  for  the  continuous  case  but  is  in  a  discrete  model,  such  as  an  array 
of  pixels.  It  represents  an  image  with  infinite  frequency  content,  meaning  that  its  Fourier 
transform  does  not  go  to  zero  as  frequency  increases  [6:20-21]. 

If  this  point  source  is  propagated  through  an  aperture,  the  psf  is  found  and  therefore  the 
transfer  function  is  also  known  for  the  system.  Figure  2.3  shows  both  the  object  and  the 
observed  image  for  several  circular  apertures.  As  the  images  show,  the  point  source 
becomes  more  and  more  spread  out  as  the  aperture  size  decreases.  This  loss  in  definition 
can  be  explained  by  the  transfer  function.  The  size  of  the  aperture  determines  the 
maximum  spatial  frequencies  that  can  be  recovered.  Since  a  point  source  contains 
infinite  frequency  content,  it  gets  blurred  out  due  to  the  aperture  cutting  off  higher  spatial 
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frequencies.  Looking  at  the  transfer  functions  in  Figure  2.4,  it  is  clear  that  as  the  aperture 
shrinks,  so  does  the  cutoff  frequency. 


1.56cm  point,  100cm  window 


Optic  =  32cm 


Optic  =  25cm  Optic  =  20cm 


Optic  =  15cm 


Optic  =  10cm 


Optic  =  5cm 


Optic  =  2cm 


Figure  2.2:  The  first  image  is  that  of  the  original  object  with  a  square  window  of  100cm2.  Each  successive  image  is  the  object 
propagated  through  a  circular  aperture  of  the  diameter  specified. 


Optic  =  20cm 


Optic  =  2cm 


Figure  2.3:  The  top  row  shows  the  image  of  a  point  source  as  seen  through  a  circular  aperture  of  the  given  diameter.  Directly  below 
each  image  is  the  transfer  function.  The  transfer  function  is  in  spatial  frequency.  As  the  diameter  shrinks,  so  does  the  maximum 
spatial  frequency. 
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2.3  Turbulence 


2.3.1  The  Long  Exposure  Optical  Transfer  Function 

The  algorithm  developed  in  this  thesis  is  based  on  the  Long  Exposure  Optical 
Transfer  Function  (LEOTF)  for  a  combined  atmosphere-optics  system.  Over  any 
integration  time,  the  optics  of  the  telescope  do  not  change  but  the  atmosphere  is  in  a 
constant  state  of  turbulent  change.  To  find  the  LEOTF  of  the  entire  system,  the 
derivation  begins  with  the  definition  of  an  OTF 


w) 


iCo  H(x,  y)H * (x  -  Afvu,  y  -  Afvv)dxdy 
y)\2dxdy 


(2.11) 


where  1  is  the  mean  wavelength,  /  is  the  effective  focal  length  of  the  optical  system,  and 
vu  and  vv  are  the  spatial  frequency  coordinates  for  the  Fourier  transfonn  space  in  the 
image  plane  [6: 139-140].  For  the  long  exposure  case  the  atmosphere  can  be  condensed 
into  a  single  phase  screen  located  directly  in  front  of  the  aperture.  With  this 
simplification 

H(x,y)  =  P(x,y)e10ix,O  (2.12) 


in  which  P(x,y )  is  the  pupil  function  and  0  is  a  random  Gaussian  variable  with  mean 
equal  to  zero  [6: 145-147].  Inserting  this  into  equation  (2. 1 1)  and  taking  the  expectation, 
the  equation  becomes 


(Vu,  ^V)  P 


P(x,y)P*(x  —  Afvu,y  —  Afvv)eld(x,y^e  xfVu-y  xfVv)dxdy 


U-Jp(x’  y)  \2dxdy 


(2.13) 


13 


where  6  is  a  random  phase  for  the  light  coming  through  the  pupil.  In  the  denominator,  the 
phase  of  the  transfer  function  is  the  same  as  the  conjugate  transfer  function  so  they  cancel 
out;  however,  this  is  not  true  for  the  numerator.  The  only  random  quantity  in  equation 
(2.13)  is  the  phase  so  the  expectation  can  be  brought  inside  upper  integrals  yielding 


-ffsys(vit<  yv) 


iCooP(x>y)p*{x  ~  tfvu> y  -  Afvv)E[eld']dxdy 

SCjp(x,y)\2dxdy 


(2.14) 


with  the  variable  substitution  O'  —  [6>(x,y)  -  0(x  -  Afvu,y  -  Afvv)]  made.  In  order  to 
detennine  a  usable  form  of  the  expectation,  the  characteristic  function  is  needed.  It  is 
defined  as 

4v(cj)  =  E[eio)Q']  (2.15) 

where  od  is  Fourier  transform  variable  associated  with  the  frequency  content  of  0'.  This 
is  similar  to  equation  (2.14)  which  has  u>  =  1.  Since  both  0(x,y)  and  6{x  -  lfvu,y  -  Ifvv ) 
are  Gaussian  random  variables,  O'  is  also  Gaussian.  The  characteristic  function  for  a  zero 
mean  Gaussian  is  given  by, 

dv(-)  =  ^  (2’16) 


where  cr2  is  the  variance  of  the  variable  0'.  Setting  oj  —  1  and  inserting  this  equivalent 
form  into  equation  (2.14)  then  pulling  it  outside  the  integrals  it  becomes 


Tfsys  (yw  yv) 


/ JCqPC*.  y)P'(x  -  Afvu,y-  Afvv)dxdy\ 
A  XCol  P(x,y)\2dxdy  ) 


(2.17) 
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a2  =  E  [{00c, y)  -  6(x  -  Afvu,y  -  A/vv)}2]  =  2 c702  -  2£’[0(x,y)0(x  -  Afvu,y  -  Afvv )]  (2.18) 


where  er2  is  the  variance  of  the  phase  difference.  Equation  (2.18)  gives  the  definition  of 
the  variance  and  then  puts  it  in  terms  of  correlations.  This  is  called  a  phase  structure 
function.  The  first  term  in  equation  (2.17)  is  only  dependent  on  the  phase  delays  caused 
by  the  atmosphere  and  the  second  term  is  just  the  OTF  of  the  optical  system  alone  which 
is  written  as 


Tf'atm  (yu>  ~^v)  ®  2 


(2.19) 


?fo(v„,  vv) 


C, p (x> y)p* (x  ~  xfvu, y  -  Afvv)dxdy 
y)\2dxdy 


(2.20) 


This  demonstrates  that  the  LEOTF  of  the  entire  system  is  completely  separable  with 
respect  to  the  atmosphere  and  optics  [5:404-407] 

-H-sys (yw  ^v)  f />;  (T; /,  Vr):}C{)(vu,  Vj,)  (2.21) 


2.3.2  Phase  Structure  Function 

One  way  to  characterize  the  statistical  correlation  between  two  points  in  a  distribution  is 
through  the  structure  function.  The  structure  function  can  be  described  as  the  variance  of 
the  difference  between  two  points  in  a  field.  For  the  case  of  the  phase  difference  the 
structure  function  looks  like 
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D(x,y,Ax,Ay)  =  E[{9(x  +  A x,y  +  Ay)  —  9(x,  y)}2] 


(2.22) 


where  x  and  y  are  general  spatial  coordinates  and  Ax  and  Ay  are  deviation  for  the 
coordinates.  At  this  point,  the  structure  function  could  be  substituted  into  equation  (2.19) 
in  place  of  the  variance  tenn.  However,  just  a  few  more  steps  will  be  performed  before 
this  is  done.  If  the  square  term  is  multiplied  out  and  the  expectation  broken  up  of  these 
tenns  it  gives 

D(x,y,  Ax,  Ay)  =  E[92(x  +  A x,y  +  Ay)]  +  E[9z(x,y)]  —  2E[9(x  +  Ax,  y  +  A y)9(x,y)]  (2.23) 

The  above  string  of  expectations  can  be  expressed  in  tenns  of  the  correlation  function 
eliminating  the  need  for  an  initial  ( x ,  y)  coordinate.  This  now  describes  the  structure 
function  as 


Dg  (Ax,  Ay)  =  2  [Rg  (0,0)  -  Re  (Ax,  Ay)]  (2.24) 

Lastly,  this  structure  function  can  be  inserted  into  equation  (2.19)  giving  an  average 
atmospheric  transfer  function  in  terms  of  autocorrelations  [1:38-40] 

Matm(Vu.Vv)  =  e~De(2X'Ay)  =  e-[K9(0,0)-Kfl(A*,A  y)]  (2.25) 

2.3.3  Kolmogorov  Structure  Function 

The  fluctuations  in  phase  are  caused  by  changes  in  the  index  of  refraction  in  the 
air.  These  changes  arise  from  fluctuations  in  temperature,  and  to  a  lesser  degree  pressure. 


16 


In  1961,  Kolmogorov  derived  an  expression  for  the  phase  structure  function.  This  is 
given  by 

D(r)  =  2.91/c2r3zC^(z)  (2-26) 

where  k  is  the  wave  number,  z  is  the  height  of  the  atmospheric  turbulence,  and  C2  is  a 
measure  of  the  strength  of  atmospheric  turbulence,  in  units  of  length.-2/3 .  Equation 
(2.26)  is  for  the  case  of  constant  C2  [5:413].  If  the  strength  is  changing  with  altitude  the 
structure  function  is  [5:428] 

D(r)  =  2.91  k2r5/3  f  C2(z')dz'  (2.27) 

”'0 


A  parameter  to  describe  atmospheric  seeing  is 


rn  d=f  0.185 


-,3/5 


[J0  c2(z')dz' \ 


(2.28) 


r0  is  called  the  Fried  Parameter  and  can  be  thought  of  as  an  average  atmospheric 
coherence  length  [5:43 1].  Rearranging  this  definition  and  replacing  the  integral  in  the 
structure  function,  it  is  written  as 

tr  \  5/3 

D(r)  =  6.88  —  (2.29) 

\r0/ 


Finally,  this  structure  function  can  be  inserted  into  the  equation  (2.25)  yielding  the  result 


-3  44(A) 

Matm(yu,vv)  =  e  VrV 


5/3 


(2.30) 
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So,  the  average  OTF  of  the  atmosphere  can  be  completely  characterized  by  a  single 
length  parameter,  r0  [5:439]. 


2.4  Detection 

The  detection  of  photons  by  a  detector  array  is  a  Poisson  process.  In  this  case,  an  event 
will  be  defined  as  an  individual  pixel  detecting  a  photon.  For  this  example,  it  will  be 
assumed  that  the  quantum  efficiency  and  photomultiplier  are  equal  to  1 .  Over  some 
period  of  time,  either  an  event  does  or  does  not  occur  with  some  probability.  A  binomial 
distribution  could  be  used  to  model  this.  However,  photon  events  occur  in  such  large 
number  that  the  number  of  trials  can  be  assumed  to  be  infinite.  In  that  case,  the 
probability  of  detection  is  given  by  a  Poisson  distribution, 


i(y)fe(y)e-i(y) 

k(y)\ 


(2.31) 


where  t(y)  is  the  mean  number  of  events  and  k(y)  is  the  number  of  successes  [5:466- 
467].  For  an  imaging  system,  the  mean  number  of  events  is  given  by  estimated  image 
intensity,  i(y),  and  k  is  the  photocount  data.  For  multiple  pixels,  the  composite  pdf  of  the 
entire  array,  assuming  statistical  independence  between  detector  measurements,  is 


G  = 


n 


j(y)ti(y)g-i(y) 

d(y)l 


(2.32) 
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where  y  is  actually  a  two  dimensional  variable  representing  the  image  plane,  Y  is  the 
dimensions  of  that  image  plane,  and  d(y)  is  the  data  [10:1065].  How  the  model 
describes  the  intensity  of  the  image  will  be  discussed  in  chapter  3. 


19 


III.  Algorithm  Design 


3.1  Log-Likelihood  Data  Model 

3.1.1  Incomplete  Data 

As  discussed  in  section  2.4,  for  photons  incident  on  a  detector,  a  Poisson  distribution  is 
used  to  model  the  random  arrival  times  of  the  light.  If  you  now  consider  the  case  of  an 
imaging  system  where  two  channels  are  statistically  independent,  the  combined  PMF  is 
just  the  multiplication  of  the  two  individual  PMFs.  This  incomplete  data  likelihood 
model  takes  the  form 


yiYi  \fiup(.yi)dupiyi)e  I,,,,(yi)^n[^(y2)dp(y2)e  lp(yz)Y 
LhA\  dUP(y  i)!  y  J 1  1 LV  dp(y2)\  ) 


where  yt  and  y2  are  general  coordinates  of  the  image  planes,  Y1  and  Y2  are  the  total 
number  of  pixels,  and  d(y)  is  the  collected  photon  counts  from  the  detectors.  For 
reasons  discussed  below,  the  data  will  be  referred  to  as  the  incomplete  data.  The 
subscript  UP  denotes  the  unpolarized  channel  and  the  subscript  P  denotes  the  polarized 
channel.  iyp(y)  and  iP(y)  are  the  average  intensity  of  the  images  given  by  the 
convolutions 

x 

tup(yi)  =  ^  0(x)h(y1  -  x)  (3.2a) 


ip(y2 )  =  ^  0(x)P(x)h(y2  -  x) 


X 


(3.2b) 
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0(x)  is  the  object,  h(y  —  x)  is  the  point  spread  function,  P(x )  is  the  set  of  polarization 
states  with  values  between  0  and  1,  and  x  is  a  general  coordinate  for  the  object  plane 
[10:1065], 


3,1.2  Complete  &  Incomplete  Data  Relationship 

At  this  point  the  concept  of  complete  data  will  be  introduced.  Complete  and  incomplete 
data  are  related  but  the  relation  between  the  two  can  be  totally  user  defined.  Further,  the 
complete  data  do  not  even  need  to  be  physical.  The  only  requirement  is  that  they  be 
statistically  consistent.  It  is  more  a  matter  of  choosing  parameters  that  make  the 
equations  solvable.  In  this  case,  the  relation  between  the  two  will  be  defined  as 

x 

d(y)  =^Td(y|x)  C3-3) 

% 


where  d(y  |x)  is  a  piece  of  complete  data  and  will  be  considered  a  Poisson  random 
variable.  This  is  saying  that  the  incomplete  data  is  the  sum  of  small  individual  pieces  of 
complete  data.  Since  the  incomplete  data  is  a  sum  over  a  set  of  Poisson  random 
variables,  it  is  also  a  Poisson  random  variable  and  remains  statistically  consistent.  The 
expectation  of  the  complete  data  is  given  by 

E[d(y\x)]  =  h(y  —  x)O(x)  (3-4) 

So,  taking  the  expectation  of  equation  (3.3)  and  inserting  equation  (3.4) 
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(3.5) 


x  x 

E[d(y)\  =  'YJE[d(y\x)]  =  X  h(y  —  x)O(x) 

X  X 


Now  the  complete  data  has  a  relationship  to  the  incomplete  data  and  is  statistically 
consistent  [10:1067]. 


3.1.3  Complete  Data  Log-Likelihood 

Taking  the  log  of  equation  (3.1)  and  assuming  statistical  independence  of  the  pixels  in  the 
images,  the  incomplete  data  log-likelihood  is 


lid='Yj  (4(y)ln  ^0(x)ft(y  —  x)j— ^0(x)h(y  -  x)  -  In  (dUP(y)\)  + 

y  \  X  '  X 

(  X  \  * 

dp(y )  In  |  0(x)P(x)fr(y  -  x)  j  -  ^  0(x)P(x)h{y  —  x)  —  In  (dP(y)!)} 


(3.6) 


The  subscript  ID  stands  for  incomplete  data.  This  is  only  shown  here  for  comparison  to 
the  complete  data  log-likelihood.  There  is  no  way  to  go  directly  from  the  incomplete  to 
the  complete  data  model.  The  complete  data  log-likelihood  is  given  by 


Lm  — 


V  Y  (dup(y)  In {0(x)h(y  -  x))  -  0{x)h{y  -  x)  -  In  (dUP(y)'.)  + 
dp{y )  ln(0(x)f’(x)/i(y  -  x))  -  0{x)P(x)h(y  -  x)  -  In  (dP(y)!)} 


(3.7) 


This  looks  very  similar  to  equation  (3.6)  but  there  is  no  direct  relationship  between  them. 
The  factorial  of  the  data  has  no  terms  for  O(x)  or  P(x)  and  will  not  affect  the  likelihood 
so  they  can  be  dropped.  Removing  these  tenns  (3.7),  it  becomes 


22 


(3.8) 


ZY  iX 

/  { dUP(y )  ln(0(x)/i(y  -  x))  -  0(x)h(y  -  x)  + 

y  'x 

dP(y)  ln(0(x)P(x)h(y  —  x))  —  0(x)P(x)h(y  —  x)} 


3.1.1  Polarization  Prior 

In  order  to  completely  describe  the  joint  likelihood  of  the  two  channel  system,  a  prior  is 
needed  to  constrain  the  polarization  parameters  P(x).  Right  now,  the  polarization  can  be 
any  number  greater  than  or  equal  to  zero,  due  to  the  Poisson  nature  of  the  signal.  It  is 
known  that  the  value  of  P  cannot  exceed  1 .  Since  the  likelihood  only  needs  to  be  limited 
by  the  upper  bound,  a  simple  “super-Gaussian”  can  be  used.  For  the  pdf,  this  would  be 

/(PO))  =  Ae~pWn,P(x )  >  0  0.9) 


where  n  is  an  even  integer  and  A  is  a  nonnalizing  factor  so  that  the  pdf  integrates  to  unity. 
The  nonnalizing  tenn  is  unimportant  since  once  the  natural  log  is  taken,  it  becomes 
another  constant  offset  like  the  data  factorial.  As  n  increases,  the  pdf  models  a  step 
function  more  and  more.  Figure  3. 1  shows  a  plot  of  this  prior  demonstrating  how 
increasing  n  approaches  a  step  function.  Its  purpose  is  to  keep  the  likelihood  low  for  P 
values  greater  than  one  but  high  for  values  less  than  one,  allowing  the  data  to  accurately 
estimate  the  P  state.  For  example,  P  =  0.95,  the  prior  would  bias  the  estimated  value 
away  from  this.  Now  the  complete  data  log-likelihood  function  is 


Lcd  =  V  V  {dUP(y)  In (O(x)ft(y  -  x))  -  0(x)ft(y  -  x)  + 
dP{y )  In (0(x)P(x)ft(y  -  x))  -  0(x)P(x)^(y  -  x)  -  P(x)n} 


(3.10) 
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Figure  3. 1 :  A  plot  of  the  prior  function  in  equation  (3.9).  As  n  increases,  the  prior  models  a  step  function  more  closely. 


3.2  Expectation  Maximization  Algorithm 

The  General  Expectation  Maximization  (GEM)  algorithm  is  a  method  for  estimating  an 
unknown  quantity,  the  complete  data,  given  a  set  of  related  known  quantities,  the 
incomplete  data.  The  process  involves  taking  the  conditional  expectation  of  the  complete 
data  likelihood  given  the  incomplete  data  and  the  old  parameter  estimates.  Once  that  is 
done,  the  final  step  involves  maximizing  the  new  function  for  the  desired  parameters. 

The  steps  for  the  GEM  algorithm  follow  the  simple  format  of: 

1.  Take  the  conditional  expectation  of  the  complete  data  log-likelihood 

2.  Maximize  the  new  equation  from  step  1  for  the  parameter(s)  of  interest 

3.  Calculate  the  conditional  expectation  given  the  old  parameters 

4.  Calculate  what  the  new  parameters  will  be  given  step  3 

5.  Iteratively  repeat  steps  3  and  4  until  a  stopping  criteria  is  reached 
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3.2.1  Expectation  Step 

Using  the  complete  data  log-likelihood  model  described  earlier  in  the  chapter,  a  new 
function  can  be  found  by  taking  the  conditional  expectation.  The  general  Q  function  is 
given  as  [10:1066] 

Q(0(x),h(y  -x))  =  E[Lcd(0(x),  h(y  —  x))|d(y)]  (3.11) 


The  Q  function  for  the  likelihood  model  in  equation  (3.10)  is  then 
Q=E[ V  V  {d„P(y)ln(0(x)/i(y  -  x))  -  0(x)h(y  -  x)  + 

'y 

dp(y)  In {0(x)P(x)h(y  —  x ))  -  0(x)P(x)h(y  —  x)  —  P(x)n }] 


From  this  point  forward,  the  conditions  of  the  expectation  will  be  suppressed  but  should 
be  kept  in  mind.  The  conditional  expectation  can  be  brought  inside  the  summation. 
However,  the  only  tenn  within  the  Q  function  that  is  a  random  variable  is  the  data. 
Making  this  simplification,  the  Q  function  becomes 


Q  =  T  y  {T[d[/p(y)  ]ln(0(x)h(y  —  x))  —  0(x)h(y  —  x)  + 

*y  ‘x 

E[dP(y)  ]ln(0(x)P(x)h(y  -  x))  -  0(x)P(x)ii(y  —  x)  —  P(x)n } 


(3.13) 


The  expectation  of  the  complete  data  is  not  a  known  quantity  and  must  be  calculated.  A 
method  for  calculating  the  expected  value  of  the  data  is  given  by  Schulz  and  used  in  the 
previous  work  done  by  Strong  [11  :Chap  5], [10: 1067].  Those  equations  are 


E[dUP(y )] 


h(y  —  x)0(x ) 

iup(y\0,h) 


dyp(y) 


(3.14a) 
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(3.14b) 


E\dP{y)] 


h(y  —  x)P(x)0(x) 
iP(y\0,P,  h ) 


dP(y) 


It  has  been  shown  that  that  the  conditional  probability  distributions  needed  to  compute 
the  conditional  expectations  in  equations  (3.14)  are  binomial  [4],  Accepting  this,  an 
intuitive  way  of  thinking  about  these  equations  is  to  remember  that  the  expectation  of  a 
binomial  distribution  is  n  times  p,  where  n  is  the  number  of  trials  and  p  is  the  probability 
of  success.  The  data  are  a  measure  of  the  number  of  photons  that  hit  a  pixel  and  can  be 
thought  of  as  n  trials.  The  numerators  in  the  equations  are  simply  a  piece  of  complete 
data.  In  the  denominator,  the  intensity  is  a  sum  over  all  pieces  of  complete  data.  So,  the 
ratio  of  these  two  terms  is  a  probability,  p. 


3.2.2  Maximization  Step 

In  the  maximization  step,  there  are  two  parameters  of  interest,  the  object,  0(x),  and  the 
polarization  state,  P(x).  These  are  solved  for  simultaneously  by  taking  the  derivative 
with  respect  to  each  and  setting  them  equal  to  zero.  First,  taking  the  derivative  of  the  Q 
function  with  respect  to  P(x0)  and  setting  equal  to  zero  gives 

y 

°  =  X  ~  -  *o)  -  nP(x o)”"1  j  (3-15) 

Note  that  the  derivative  was  taken  with  respect  to  a  specific  instance  of  x  so  the 
summation  over  the  object  plane  vanishes.  The  object  in  the  second  tenn  is  not 
dependent  on  y,  only  on  the  psf.  By  definition,  assuming  and  ideal  optical  system,  the 
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sum  of  the  psf  is  equal  to  1.  Including  this  and  solving  for  O(x0)  the  maximization 
becomes 


0  Oo) 


^g[dp(y)] 

P(x. o) 


^  nP(x0)?l  1 


(3.16) 


This  is  as  far  as  this  equation  needs  to  be  simplified.  Now  to  repeat  this  procedure  with 
respect  to  the  object,  O(x0).  Taking  the  derivative  and  setting  equal  to  zero,  the 
maximization  is 


n^Vf^W] 
4n  °(*o) 


-  h(y  -  x0)  + 


E[dP(y)] 
O(x0 ) 


P(x0)h(y  -  x0 ) 


(3.17) 


Again,  summing  over  the  psf  terms  and  solving  for  O(x0),  the  final  form  is 


O(x0)  = 


ZyEjdupiy)]  +T^E[dP(y')] 
1  +  P(x0) 


(3.18) 


To  remove  O(x0),  equations  (3.16)  and  (3.18)  are  set  equal  to  each  other, 


JZE[dUP(y)]  +  'ZyE[dP(y)]  _'ZlE[dp(y)]  ^  ^  ^n_1 

iTp(^  "  — po£ j  Z nP(Xo) 


(3.19) 


Moving  all  tenns  to  one  side  and  separating  them  gives 


'ZyE[dUP(y)] 
1  +  P(x0) 


g p[[dp00] 
1  +  P(x0) 


y  y 

dp(y)  -^nP(x0)n_1 

y  y 


(3.20) 
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The  last  step  is  to  multiply  through  by  —  (1  +  P(x0))  which  shows  that 

Y  Y 

2/[dP(y)]  (3.21) 

y 

At  this  point,  the  only  way  to  solve  for  P  (x)  is  to  find  the  roots  of  this  polynomial 
expression.  As  mentioned  before,  as  n  grows  larger,  the  more  accurately  the  prior  models 
a  step  function.  The  drawback  to  a  large  n  is  that  it  requires  more  computing  time.  There 
is  only  one  real  root  to  the  expression,  but  which  one  it  is  can’t  be  known  prior  to  solving 
for  them  all.  A  practical  value  to  choose  for  n  would  be  a  number  less  than  10. 

In  the  algorithm,  once  P(x)  is  estimated,  it  is  inserted  back  into  equation  (3.18).  Then, 
an  estimate  for  the  object  can  be  found.  Having  a  new  object  and  polarization  state,  the 
algorithm  is  looped  so  that  a  new  intensity  and  complete  data  expectation  can  be 
calculated  to  find  the  next  P(x)  and  0(x)  estimates.  The  iteration  process  is  repeated 
until  the  stopping  criterion  is  reached. 


0  =  nP(x)n+1  +  nP(x)”  + 


P(x)^E[dUP 


(y)] 


3.2.3  Stopping  Criterion 

In  order  to  know  when  the  algorithm  should  be  stopped,  a  specific  criterion  needs  to  be 
set.  One  way  to  do  this  is  to  use  the  statistical  properties  of  the  noise.  The  intensity  is  a 
convolution  of  the  estimated  object  with  the  total  system  transfer  function;  in  effect,  it 
tries  to  mimic  the  data.  However,  the  recreated  object  cannot  exactly  account  for  the 
random  noise  in  the  system.  So,  once  the  variance  difference  between  the  intensity  image 
and  the  data  comes  within  one  standard  deviation  of  the  data’s  variance,  the  algorithm  is 


28 


stopped.  If  it  were  allowed  to  iterate  beyond  this  point,  the  reconstruction  will  try  to 
incorporate  too  much  of  the  noise  as  part  of  the  object,  skewing  the  object  estimate.  The 
stopping  criterion  is  then  [8] 

Y 

Op  <  2jdP(y)  -  iP(y)]2  (3-22) 

y 

The  polarized  data  set  is  used  instead  of  the  unpolarized  because  the  variance  is  of  the 
polarized  image  is  always  less  and  the  algorithm  will  converge  faster.  In  test  runs  of  the 
algorithm,  when  the  unpolarized  data  was  used,  the  estimated  object  was  noticeably  over¬ 
iterated. 


3.3  Seeing  Parameter  Estimation 

The  true  seeing  parameter  r0  value,  which  is  used  to  generate  the  psf,  is  not  known  when 
the  algorithm  is  initialized.  The  reconstruction  is  simply  recalculated  for  a  range  of  ro 
values.  The  likelihood  of  the  reconstructed  image  is  found  each  time  by  using  equation 
(3.8).  After  processing  is  done,  a  likelihood  curve  vs.  ro  values  can  be  plotted.  In  all 
cases  tested,  this  plot  looks  similar  to  Figure  3.2.  The  likelihood  sharply  rises  from  very 
bad  seeing  to  the  correct  ro  value.  This  gives  a  maximum  near  the  actual  seeing 
parameter.  The  likelihood  value  itself  is  not  important;  it  is  simply  a  comparative 
measure  between  different  ro  values.  The  algorithm  has  a  tendency  to  slightly 
overestimate  ro  but  is  usually  within  a  half  centimeter  of  truth. 
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Likelihood 


Figure  3.2:  After  the  algorithm  produces  an  estimated  image  for  a  range  of  revalues,  the  likelihood  is  calculated  and  plotted.  It 
reaches  maximum  at  or  near  the  true  seeing  parameter.  In  this  case,  the  revalue  was  6cm  and  was  estimated  exactly. 
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IV.  Results 


4.1  Simulation  Setup 

In  order  to  accurately  quantify  the  capabilities  of  the  polarization  diversity  algorithm,  a 
simple  object  is  needed  to  be  reconstructed.  It  is  not  as  impressive  as  using  satellite 
images  but  will  provide  substantially  more  infonnation  about  the  bounds  of  the 
algorithm.  Two  bar  targets  separated  by  two  pixels  are  used  as  the  object.  The  total 
image  size  used  is  60cm  on  a  64x64  pixel  grid  so  each  pixel  represents  0.94cm.  This 
object  is  shown  if  Figure  4. 1 .  Each  bar  is  set  with  a  different  polarization  state.  The 
following  sections  describe  the  tests  that  were  done  on  the  bar  target  by  varying  different 
parameters  of  the  data. 


Bar  Target  Object 
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30 1 

Pixels 

40 1 


50 


6°yyrrp4  h| 

10  20  30  40  50  60 

Pixels 

Figure  4. 1 :  Two  bars  to  be  propagated  through  a  simulated  atmosphere  and  telescope  aperture.  Each  bar  is  given  a  different 
polarization  parameter. 
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4.2  Test  One:  Signal  Variation 

The  first  simulations  were  done  by  varying  the  polarizations  on  both  bars  and  also 
increasing  the  signal  strength  of  the  data.  P  for  both  bars  was  iterated  from  0  to  1  in  0.1 
steps  separately.  Once  there  was  an  object  estimate  for  all  possible  polarization 
combinations,  the  signal  power  was  increased  and  the  process  was  repeated. 

To  keep  data  sizes  small,  a  one-dimensional  slice  is  taken  for  each  estimated  object,  since 
the  object  is  uniform  in  the  vertical  dimension.  If  the  object  were  reconstructed  perfectly, 
it  would  look  like  two  periods  of  a  square  wave,  as  in  Figure  4.2.  In  order  to  do  that 
though,  infinite  frequency  content  would  be  required  in  the  continuous  case. 
Unfortunately  too  much  of  the  higher  spatial  frequencies  are  lost  through  diffraction. 
However,  the  goal  of  this  simulation  is  not  to  exactly  reproduce  the  bar  targets  but  rather 
just  to  distinguish  that  there  are  two  targets  present.  An  example  of  the  reconstructed 
object  is  shown  in  Figure  4.3.  It  should  be  noted  that  the  intensity  of  the  image  is  much 
higher  than  that  of  the  object.  This  does  not  matter  as  they  are  both  scaled  to  the  same 
intensity  range  when  displayed.  Figure  4.4  is  a  screenshot  of  the  algorithm  in  progress 
for  bar  polarizations  of  0.2  and  1.0.  The  top  two  images  are  that  of  the  observed  data 
channels  after  being  propagated  through  the  atmosphere  and  aperture.  The  bottom  right 
is  the  object  and  the  bottom  left  is  the  reconstruction. 
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Figure  4.4:  A  screenshot  of  the  algorithm  reconstruction  the  object.  The  upper  left  shows  the  unpolarized  data  channel,  polarized 
channel  in  the  upper  right,  object  estimate  in  lower  left,  and  object  in  lower  right. 


To  characterize  how  well  the  individual  bars  can  be  resolved,  a  ratio  is  taken  of  the 
minimum  between  the  two  bars  and  the  maximum  of  the  bars.  This  gives  a  number  less 
than  one  with  smaller  ratios  indicating  better  resolution.  To  see  how  different 
polarizations  affect  resolution,  a  contour  plot  of  the  ratios  is  shown  in  Figure  4.5.  The  x 
and  v  axis  are  the  polarizations  of  the  two  bars.  The  figure  shows  that  the  best  resolutions 
are  found  where  at  least  one  of  the  targets  is  highly  polarized,  close  to  either  0  or  1 .  The 
area  near  zero  polarization  diversity,  along  the  diagonal,  does  not  fare  as  well.  Along 
with  finding  where  the  best  resolutions  occur,  the  average  value  can  be  found  with 
respect  to  the  signal  level. 
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Polarization  Contour  Plot 


Figure  4.5:  An  example  contour  plot  of  resolving  capability  based  on  polarization  of  bar  targets 


Test  2:  Reproducibility 

The  next  simulation  produced  the  same  data  as  test  one  but  ran  only  two  signal  levels, 

2x  104  and  9X 104  photons,  with  an  rg  of  3.0cm.  The  aperture  size  for  all  tests  was  set  to 
10cm.  However,  this  was  repeated  5  times  each  to  make  sure  that  similar,  not  necessarily 
identical,  contour  plots  were  produced.  The  5  contour  plots  were  then  averaged.  Most  of 
the  center  areas  where  there  is  not  much  diversity  do  not  show  comparatively  high 
resolutions.  However,  near  the  very  edges  where  diversity  is  the  highest,  the  resolutions 
are  much  better.  This  showed  that  similar  results  were  achieved  each  time  the  simulation 
was  initiated. 
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Test  3:  Resolvable  Threshold 


Through  testing  a  set  of  different  ro  values,  aperture  size  of  10cm,  the  threshold  at  which 
the  algorithm  can  distinguish  two  different  objects  was  found  to  be  2.2cm.  At  this  value, 
high  polarization  diversity  yields  some  resolution  but  not  much  whereas  low  diversity 
still  cannot  resolve  the  two  objects.  If  ro  is  incremented  to  2.3cm,  resolution  increases 
dramatically  across  all  polarization  combinations.  Figure  4.6  shows  a  very  nice 
symmetric  gradient  across  all  the  plots.  As  the  diversity  decreases  toward  the  center  line, 
resolution  gets  worse  and  then  begins  to  get  better  as  it  approaches  the  other  extreme.  It 
should  be  fairly  symmetric  since  which  bar  is  polarized  a  certain  way  shouldn’t  affect  the 
reconstruction. 

The  first  run  of  this  test  was  done  at  a  5  x  1 04  photocount.  Repeating  the  same  simulation 
with  lower  signal  values  of  10  ,5x10,10,  and  500  (r0  of  2.3cm),  yields  expectedly 
worse  results  along  the  diagonal  where  there  is  no  diversity.  In  contrast,  as  the  diversity 
increases,  the  signal  level  becomes  less  and  less  important  and  the  reconstruction 
continues  to  produce  very  good  results.  At  a  certain  point  below  a  photocount  of  100,  the 
method  for  quantifying  the  data  breaks  down.  The  reconstructed  image  is  no  longer 
symmetric.  The  cross  section  does  not  accurately  represent  the  image  as  a  whole. 
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Figure  4.6:  Contour  plot  of  resolving  capability  based  on  bar  polarizations.  Signal  levels  shown  are  500  (top  left),  1000  (top  center), 
5000  (top  right),  10000  (bottom  left),  50000  (bottom  right) 


Test  4:  Sample  Satellite  Reconstruction 

This  algorithm  was  designed  to  reconstruct  images  of  satellites  acquired  from  ground- 
based  telescopes.  In  order  to  make  a  qualitative  estimate  of  how  well  the  algorithm 
works,  a  simulated  satellite  is  used.  Figure  4.7  shows  the  unpolarized  satellite  and  the 
polarized  satellite  before  propagation.  The  polarized  data  was  created  using  an  edge 
detection  filter  on  the  unpolarized  object  to  represent  the  high  polarization  of  edges.  This 
is  not  the  best  representation  of  a  polarized  satellite  but  it  is  sufficient  for  this  test.  Figure 
4.8  shows  another  screenshot  of  the  algorithm  with  both  data  channels,  the  object  and  the 
reconstruction.  There  is  noticeably  better  resolution  on  the  reconstruction  than  the  raw 
data. 
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Figure  4.7:  A  satellite  with  the  unpolarized  object  (left)  and  polarized  object  (right) 


Figure  4.8:  Screenshot  of  the  algorithm  restoring  a  satellite  image.  Unpolarized  data  channel  (top  left),  polarized  data  channel  (top 
right),  object  estimate  (bottom  left),  object  (bottom  right) 
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V.  Conclusions 


5.1  General  Conclusions 

The  purpose  of  the  algorithm  developed  in  this  thesis  was  to  produce  improved  object 
estimates  beyond  that  of  single  channel  models  currently  employed.  The  process  is 
dependent  on  not  just  polarization  but  polarization  diversity.  When  two  points  in  the 
image  have  significantly  different  polarizations,  one  horizontal  and  one  vertical  for 
instance,  the  ability  of  the  algorithm  to  resolve  them  is  greatly  increased.  In  contrast  to 
this,  when  there  is  little  to  no  diversity  the  algorithm  will  only  resolve  as  well  as  a  single 
channel  system  would. 

The  first  unexpected  consequence  of  using  polarization  diversity  is  that  it  is  invariant  to 
signal  level.  When  there  is  high  diversity,  low  signal  levels  can  be  resolved  just  as  well 
as  much  higher  ones.  This  continues  down  to  SNR’s  below  10.  However,  as  the  signal 
decreases,  points  with  low  polarization  diversity  cannot  be  distinguished  as  well  or  at  all. 
The  required  level  of  diversity  increases  as  the  signal  fades.  The  points  must  be 
increasingly  polarized  to  achieve  the  same  resolution  but,  at  the  extremes  points,  they  are 
always  resolved  to  the  same  high  level. 

The  other  goal  of  the  algorithm  was  to  determine  the  Fried’s  parameter  for  the 
atmosphere  without  having  any  knowledge  of  the  seeing  conditions.  The  work  done  by 
MacDonald  [7:5-(29-30)j  required  a  secondary  method  to  calculate  the  ro  value  from  the 
likelihood  curve.  Without  a  separate  estimation  of  the  ro  value,  the  likelihood  curve  rises 
steeply  and  levels  off  but  never  reaches  an  apex.  In  the  case  of  polarization  diversity,  this 
additional  calculation  is  not  needed.  The  curve  does  reach  a  maximum  and  then  begins  to 
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decrease.  With  the  additional  information  provided  by  a  two  channel  system,  the  correct 
seeing  parameter  naturally  falls  out  of  the  likelihood. 


5.2  Future  Work 

One  simulation  that  was  not  preformed  was  to  detennine  the  exact  resolving  capability 
based  on  polarization  and  bar  separation.  This  simulation  would  show  how  well  or  badly 
the  algorithm  can  resolve  compared  to  the  diffraction  limit  of  the  optic  used.  Previous 
work  done  by  Strong  showed  that  resolutions  up  to  twice  the  diffraction  limit  are  possible 
using  the  polarization  diversity  method  he  developed  [11]. 

The  system  used  in  this  thesis  is  two-channel,  with  one  unpolarized  and  one  linearly 
polarized.  The  algorithm  could  be  expanded  to  be  a  multi-channel  system  incorporating 
several  polarization  parameters,  including  the  other  linear  polarization  and  both  left  and 
right  circularly  polarized  states.  Adding  in  the  third  channel  for  the  other  linear  state  is 
relatively  easy.  The  circular  polarizations  would  require  a  greater  effort. 

The  algorithm  has  yet  to  be  tested  using  measured  data.  A  suggested  method  for  a 
controlled  test  is  to  use  a  small  telescope  pointed  towards  a  bar  target  along  a  horizontal 
path.  Just  as  a  simulated  bar  target  was  used  in  chapter  4,  a  physical  one  could  be  used 
with  a  polarizer  placed  in  front  of  each  one  of  the  bars.  Using  a  horizontal  path  can  allow 
for  the  average  ro  value  to  be  accurately  estimated.  Images  could  then  be  taken  through 
the  telescope  and  processed  using  the  algorithm  in  the  exact  same  way  that  the  simulated 
data  was. 
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Two  observatory  sites  are  currently  attempting  to  resolve  Geo  satellites,  University  of 
Arizona  MMT,  Figure  5.1,  and  the  Maui  Space  Surveillance  Center,  Figure  5.2.  The 
algorithm  could  be  tested  and  implemented  for  data  acquired  at  these  sites.  With  the 
notable  degree  of  polarization  usually  found  in  reflected  light  from  satellites  and  the  low 
signal  levels  from  them,  it  is  a  prime  candidate  for  this  algorithm. 


Another  possible  application  is  in  Solar  Physics.  The  Sun’s  intense  magnetic  field  has  an 
effect  on  the  light  emitted  from  the  solar  surface.  The  Zeeman  Effect  splits  light  emitted 
in  strong  magnetic  fields.  A  slight  shift  in  wavelength  occurs  slightly  above  and  below 
the  non-degenerate  state.  The  separate  spectroscopic  lines  are  also  polarized.  This 
polarization  is  very  noticeable  especially  in  and  near  sunspot  regions.  The  polarization 
diversity  algorithm  can  be  used  to  resolve  these  sunspot  features.  Due  to  the  high  signal 
levels  received  from  the  sun,  adaptive  optics  are  not  usually  feasible  in  solar  observation 
making  this  an  ideal  application  [2:150-152,418-419]. 


Figure  5.1:  6.5  meter  primary  telescope  mirror  at  the 
University  of  Arizona  MMT. 


Figure  5.2:  Air  Force  Maui  Optical  Station.  Haleakala,  Maui,  Hawaii . 

Image  take  from  AMOS  website: 

http ://  www.  maui.  afmc .  af.mil/Photos/MS  S  S2_lg.j  pg 
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